#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <pthread.h>

#define PREC 35000000

mpf_t aBef, aAft, bBef, bAft, tBef, tAft, pBef, pAft, pi, aux, aux_a, aux_b, aux_p, aux1, aux2, aux3;

void *calc_aAft(void *x)
{
   mpf_add(aux_a,aBef,bBef);
   mpf_div_ui(aAft,aux_a,2);   // aAft = (aBef + b)/2;
   pthread_exit(NULL);
}

void *calc_b(void *x)
{
   mpf_mul(aux_b,aBef,bBef);
   mpf_sqrt(bAft,aux_b);        // b =  sqrt(aBef*b);
   pthread_exit(NULL);
}

void *calc_t(void *x)
{
   mpf_sub(aux2,aBef,aAft);
   mpf_pow_ui(aux1,aux2,2);
   mpf_mul(aux2,pBef,aux1);
   mpf_sub(aux3,tBef,aux2);  	   //t = t - p*pow((aBef - aAft),2);
   mpf_set(tAft,aux3);
   pthread_exit(NULL);
}

void *calc_p(void *x)
{
   mpf_mul_ui(aux_p,pBef, 2);    //p = 2*p;
   mpf_set(pAft, aux_p);
   pthread_exit(NULL);
}

int main(void)
{
   int i, ret;
   pthread_t thread_aAft, thread_b, thread_t, thread_p;
   pthread_attr_t attr;

   pthread_attr_init(&attr);
   pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

   mpf_init2(aBef, PREC);
   mpf_init2(aAft, PREC);
   mpf_init2(bBef, PREC);
   mpf_init2(bAft, PREC);
   mpf_init2(tBef, PREC);
   mpf_init2(tAft, PREC);
   mpf_init2(pBef, PREC);
   mpf_init2(pAft, PREC);
   mpf_init2(pi, PREC);
   mpf_init2(aux, PREC);
   mpf_init2(aux_a, PREC);
   mpf_init2(aux_b, PREC);
   mpf_init2(aux_p, PREC);
   mpf_init2(aux1, PREC);
   mpf_init2(aux2, PREC);
   mpf_init2(aux3, PREC);

   mpf_set_ui(pBef,1);
   mpf_set_ui(aBef,1);

   mpf_sqrt_ui(aux,2);
   mpf_ui_div(bBef,1,aux);  //iniciando b

   mpf_set_ui(aux,1);	
   mpf_div_ui(tBef,aux,4);    //iniciando t

   for(i = 0;i < 24;i++)
   {
      ret = pthread_create(&thread_aAft, &attr, calc_aAft, NULL);
      if(ret) exit(-1);
      ret = pthread_create(&thread_b, &attr, calc_b, NULL);
      if(ret) exit(-1);
      ret = pthread_create(&thread_p, &attr, calc_p, NULL);
      if(ret) exit(-1);
      ret = pthread_join(thread_aAft, NULL);
      if(ret) exit(-1);
      ret = pthread_create(&thread_t, &attr, calc_t, NULL);
      if(ret) exit(-1);
      ret = pthread_join(thread_b, NULL);
      if(ret) exit(-1);
      ret = pthread_join(thread_t, NULL);
      if(ret) exit(-1);
      ret = pthread_join(thread_p, NULL);
      if(ret) exit(-1);
      mpf_set(aBef,aAft);	  //aBef = aAft;
      mpf_set(bBef,bAft);	  //bBef = bAft;
      mpf_set(tBef,tAft);	  //tBef = tAft;
      mpf_set(pBef,pAft);	  //pBef = pAft;
   }

   mpf_add(aux2,aBef,bBef);
   mpf_pow_ui(aux1,aux2,2);
   mpf_set(aux3,tBef);
   mpf_mul_ui(tBef,aux3,4 );
   mpf_div(pi,aux1,tBef);  	 //pi = (pow((aBef + b),2))/(4*t);

   gmp_printf("Pi = %.100Ff\n",pi);

   mpf_clear(aux);
   mpf_clear(aux_a);
   mpf_clear(aux_b);
   mpf_clear(aux_p);
   mpf_clear(aux1);
   mpf_clear(aux2);
   mpf_clear(aux3);
   mpf_clear(aBef);
   mpf_clear(aAft);
   mpf_clear(bBef);
   mpf_clear(tBef);
   mpf_clear(pBef);
   mpf_clear(bAft);
   mpf_clear(tAft);
   mpf_clear(pAft);
   mpf_clear(pi);

   pthread_exit(NULL);
}
